home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_4.3 / xah / xah.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  12.0 KB  |  522 lines
  1. /* 
  2.  * xah.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * X(Extract) A(rea) H(istogram)
  11.  *
  12.  */
  13. #include "xah.h"
  14.  
  15. #define    ON            1
  16. #define    OFF            0
  17. #define    DISABLE            -1
  18.  
  19. #define    PROMPT_CUTOFF        OFF      /* prompt for cutoff value(ON/OFF) */
  20. #define    DISPL_AOI        OFF
  21. #define    SCL_AOI            OFF
  22. #define    BIN_AOI            ON            /* binarize separate AOI (ON) or */
  23.       /* binarize and scan same AOI (OFF) */
  24.       /* disable binarization */
  25. #define    LIM_AREA        ON            /* limit range of dom areas */
  26. #define    WRITE_FILE        ON          /* prompt for outp file name */
  27. #define    WRITE_ALL        ON           /* dump bub data to addit outp file */
  28. #define    IP_TOOL            ON
  29.  
  30.  
  31. #define MAX_ND            10000          /* max permiss. number of domains */
  32. #define    NSQ             6                /* histogram limited to NSQ*sigma */
  33.  
  34.  
  35. #define    DEFAULT            0             /* scale hist by min and max data */
  36. #define    FIX_BINWD        1
  37. #define    FIX_RANGE        2            /* restr. range, center at mean */
  38. #define    FIX_INTERVAL        3         /* preset min AND max data */
  39.  
  40. #define    MIN_ND            5              /* min no domains for area histog   */
  41.  
  42. #define    BUF_SZ            256            /* sz of input file buffer */
  43.  
  44.  
  45. #define    SHOW_HIST
  46. #define    DBG_HIST
  47. #define    DEBUG
  48. #define    S_DEBUG
  49.  
  50.  
  51. /*
  52.  *  Global Variables
  53.  */
  54. static char HffPath[129];       /* path for HFF files            */
  55. static char Scratch[32];        /* Scratch pad temporary buffer  */
  56.  
  57. static char *fn_prompt =
  58. {"...filename (.hdt)?"};
  59. static char hbuf[BUF_SZ], *phbuf = &hbuf[0];  /* file name buf for hdt output file */
  60. static char vbuf[BUF_SZ], *pvbuf = &vbuf[0];  /* file name buf for vin output file */
  61. static char gbuf[BUF_SZ], *pgbuf = &gbuf[0];  /* file name buf for gdt output file */
  62. extern char *optarg;
  63. extern int optind, opterr;
  64. extern short tiffInput;         /* flag=0 if no ImageIn to set tags; else =1 */
  65.  
  66.  
  67. int SORT_STAT = 1;              /* flag to ind sorted inp to eval_hist() */
  68.  
  69. int CDYN = 1;                   /* flags ind data type */
  70. int X_SQ = 1;                   /* when set, eval hist of SQRT(area/mean_a) */
  71.  
  72.  
  73. /*
  74.  * error message
  75.  */
  76. void
  77. exitmess (prompt, status)
  78.      char *prompt;
  79.      int status;
  80. {
  81.   printf (prompt);
  82.   printf ("\n");
  83.  
  84.   exit (status);
  85. }
  86.  
  87.  
  88. void
  89. fail_alloc (char *str, int code)
  90. {
  91.   printf ("\n...memory alloc for %s failed\n", str);
  92.   exit (code);
  93. }
  94.  
  95.  
  96.  
  97. /*
  98.  * comparison of two elements in area array (for qsort)
  99.  * here: array elements are of type unsigned int
  100.  */
  101. int
  102. compare (a1, a2)
  103.      int *a1;
  104.      int *a2;
  105. {
  106.   return ((int) SIGN (*a1 - *a2));
  107. }
  108.  
  109.  
  110. /* sort bubble structure on y, then x, coord */
  111. int
  112. bcomp (b1, b2)
  113.      struct bubble *b1;
  114.      struct bubble *b2;
  115. {
  116.  
  117.   if (b1->ctr.y < b2->ctr.y)
  118.     return (-1);
  119.   if (b1->ctr.y > b2->ctr.y)
  120.     return (1);
  121.   if (b1->ctr.x < b2->ctr.x)
  122.     return (-1);
  123.   if (b1->ctr.x > b2->ctr.x)
  124.     return (1);
  125.  
  126.   return (0);
  127. }
  128.  
  129.  
  130.  
  131. /*
  132.  * usage of routine
  133.  */
  134. void
  135. usage (char *progname)
  136. {
  137.   progname = last_bs (progname);
  138.   printf ("USAGE: %s inimg outimg [-a x1 y1 x2 y2] [-L]\n", progname);
  139.   printf ("\n%s constructs area moments and area histogram,\n", progname);
  140.   printf ("of image containing domains (\"blobs\");\n");
  141.   printf ("the output image is a point pattern representing\n");
  142.   printf ("the domain centroids; also evaluates histogram of\n");
  143.   printf ("domain areas\n\n");
  144.   printf ("ARGUMENTS:\n");
  145.   printf ("       inimg: input image filename (TIF)\n");
  146.   printf ("      outimg: output image filename (TIF)\n\n");
  147.   printf ("OPTIONS:\n");
  148.   printf ("   -a x1 y1 x2 y2: upper left(x1, y1), lower right(x2, y2)\n");
  149.   printf ("                   coords of area to scan (int);\n");
  150.   printf ("                   default is area of imgin\n");
  151.   printf ("               -L: print Software License for this module\n");
  152.   exit (1);
  153. }
  154.  
  155. void
  156. main (int argc, char **argv)
  157. {
  158.   struct bubble *ocbub, *rcbub, *cbub;
  159.   unsigned int *oarea, *rarea, *area;
  160.  
  161.   double mean_a, sdev_a, skew_a;
  162.   double del_a;
  163.   int id, nd, ird;
  164.   int c;
  165.  
  166.   float *farea;
  167.  
  168.  
  169.   int BIN_METHOD = DEFAULT;     /* select binning criterion */
  170. /*      float *ahst; */
  171.   struct Histo *ah;
  172.  
  173.   double foo = 1.0;
  174.   int nparms = 3;
  175.   FILE *wfp, *vfp, *gfp;
  176.   Image *imgIn;
  177.   Image *imgOut;
  178.   int aoi_height;
  179.   int aoi_width;
  180.   int i_arg;
  181.   int x1 = -1, y1 = -1, x2 = -1, y2 = -1;
  182.   char in_buf[IN_BUF_LEN];
  183.  
  184. /*
  185.  * cmd line options
  186.  */
  187.   static char *optstring = "a:L";
  188.  
  189. /*
  190.  * parse command line
  191.  */
  192.   optind = 3;                   /* set getopt to point to the start of the arg list */
  193.   opterr = ON;                  /* give error messages */
  194.  
  195.  
  196.   if (argc < 3)
  197.     usage (argv[0]);
  198.  
  199.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  200.     switch (i_arg) {
  201.     case 'a':
  202.       if (!sscanf (argv[optind - 1], "%d", &x1) || !sscanf (argv[optind], "%d", &y1)
  203.           || !sscanf (argv[optind + 1], "%d", &x2) || !sscanf (argv[optind + 2], "%d", &y2)) {
  204.         printf ("Error getting values for area to scan\n");
  205.         printf ("Will attempt to scan entire picture\n");
  206.       }
  207.       break;
  208.     case 'L':
  209.       print_sos_lic ();
  210.       exit (0);
  211.     default:
  212.       printf ("\ngetopt: unknown condition encountered\n");
  213.       exit (1);
  214.       break;
  215.     }
  216.   }
  217.  
  218. /*
  219.  * Read input image
  220.  */
  221.   imgIn = ImageIn (argv[1]);
  222.  
  223.   if (imgIn->bps == 8 && imgIn->spp == 3) {
  224.     printf ("Got RGB image!!!\nInput image must be Grayscale or B&W!!\n");
  225.     exit (1);
  226.   }
  227.   if (x1 == -1 || y1 == -1 || x2 == -1 || y2 == -1) {
  228.     x1 = 0;
  229.     y1 = 0;
  230.     x2 = imgIn->width - 1;
  231.     y2 = imgIn->height - 1;
  232.   }
  233.   else if (x1 < 0 || x1 >= imgIn->width
  234.            || y1 < 0 || y1 >= imgIn->height
  235.            || x2 < 0 || x2 >= imgIn->width
  236.            || y2 < 0 || y2 >= imgIn->height
  237.            || x1 > x2 || y1 > y2) {
  238.     printf ("Invalid values for area to scan\n");
  239.     printf ("Setting to size of input image\n");
  240.     x1 = 0;
  241.     y1 = 0;
  242.     x2 = imgIn->width - 1;
  243.     y2 = imgIn->height - 1;
  244.   }
  245.   aoi_height = y2 - y1 + 1;
  246.   aoi_width = x2 - x1 + 1;
  247.  
  248. /*
  249.  * Allocate memory for output image
  250.  */
  251.   imgOut = ImageAlloc (aoi_height, aoi_width, imgIn->bps);
  252.  
  253.  
  254.   /*
  255.    * Draw a 255 border around image to eliminate
  256.    * edge effects for line fills
  257.    */
  258.   draw_rect (0, 0, imgOut->width - 1, imgOut->height - 1, imgIn, WHITE);
  259.  
  260.   /* 
  261.    * measure and record domain areas, mark and record domain centroids
  262.    */
  263.   if ((ocbub = (struct bubble *) calloc (MAX_ND, sizeof (struct bubble))) == NULL)
  264.       fail_alloc ("ocbub", 1);
  265.   nd = ah_AOI (x1, y1, x2, y2, imgIn, imgOut, ocbub, MAX_ND);
  266.  
  267.   ImageOut (argv[2], imgOut);
  268.   if ((oarea = (unsigned int *) calloc (MAX_ND, sizeof (unsigned int))) == NULL)
  269.       fail_alloc ("oarea", 1);
  270.  
  271. /*
  272.  * sort area array
  273.  */
  274.  
  275.   if (nd > MIN_ND) {
  276.     extract_area (nd, oarea, &ocbub[0]);
  277. #ifdef DEBUG
  278.     for (id = 0; id < nd; id++) {
  279.       printf ("  %6u", (ocbub + id)->area);
  280.       if ((id + 1) % 8 == 0)
  281.         printf ("\n");
  282.     }
  283. #endif
  284.  
  285.   }
  286.   else {
  287.     printf ("\n...number of domains below lower limit of %d\n",
  288.             MIN_ND);
  289.  
  290.     exit (1);
  291.   }
  292.  
  293.  
  294. /*
  295.  * evaluate moments
  296.  */
  297.  
  298. #ifdef DBG_HIST
  299.   printf ("\n...evaluating mean area (area array):");
  300. #endif
  301.  
  302.   mean_a = eval_mean (nd, oarea);
  303.   printf (" (raw) mean area: %6.2lf", mean_a);
  304.  
  305.  
  306. #ifdef DBG_HIST
  307.   printf ("\n...evaluating rmsq (area array):");
  308. #endif
  309.   sdev_a = eval_sdev (mean_a, nd, oarea);
  310.   printf ("...rmsq dev: %6.2lf", sdev_a);
  311.  
  312.  
  313. #ifdef DBG_HIST
  314.   printf ("\n...evaluating skew (area array):");
  315. #endif
  316.   skew_a = 0.0;
  317.   if (sdev_a > 0.01) {
  318.     skew_a = eval_skew (mean_a, sdev_a, nd, oarea);
  319.     printf ("...skew: %6.2lf", skew_a);
  320.   }
  321.   printf ("\n");
  322.  
  323.  
  324.  
  325. /*
  326.  * in imgOut, generate a display of the pattern marking sites
  327.  * by UPW_TRIANGLE if the corresponding area exceeds mean_a, and
  328.  * by a DNW_TRIANGLE otherwise
  329.  */
  330.   encode_ba (mean_a, sdev_a, ocbub, nd, imgOut, GRAY);
  331.  
  332. /* 
  333.  * limit values in array oarea to those within a spread of NSQ*sdev_a
  334.  * to eliminate artefacts
  335.  */
  336.   area = oarea;
  337.   cbub = ocbub;
  338.   if (LIM_AREA == ON) {
  339.     if ((rarea = (unsigned int *) calloc (nd, sizeof (unsigned int))) == NULL)
  340.         fail_alloc ("rarea", 1);
  341.  
  342.     if ((rcbub = (struct bubble *) calloc (nd, sizeof (struct bubble))) == NULL)
  343.         fail_alloc ("rcbub", 1);
  344.  
  345.  
  346.     ird = 0;
  347.     for (id = 0; id < nd; id++) {
  348.       del_a = (double) (*(oarea + id)) - mean_a;
  349.       if (fabs (del_a) <= NSQ * sdev_a) {
  350.         *(rarea + ird) = *(oarea + id);
  351.         *(rcbub + ird) = *(ocbub + id);
  352.         ird++;
  353.       }
  354.     }
  355.     nd = ird;
  356.     printf ("...%d domains within %d*std dev of mean area:\n",
  357.             nd, NSQ);
  358.  
  359. /*
  360.  * handle memory alloc
  361.  */
  362.     rarea = (unsigned int *) realloc (rarea, nd * sizeof (unsigned int));
  363.     rcbub = (struct bubble *) realloc (rcbub, nd * sizeof (struct bubble));
  364.  
  365.     area = rarea;
  366.     cbub = rcbub;
  367.  
  368.  
  369. #ifdef DEBUG
  370.     printf ("...domain areas within %d*rmsq of mean...\n", NSQ);
  371.     for (id = 0; id < nd; id++) {
  372.       printf ("  %6u", (cbub + id)->area);
  373.       if ((id + 1) % 8 == 0)
  374.         printf ("\n");
  375.     }
  376.     printf ("\n");
  377. #endif
  378.   }
  379. /*
  380.  * bin data into histogram
  381.  */
  382. #ifdef DBG_HIST
  383.   printf ("\n...evaluate area histogram...\n");
  384. #endif
  385.  
  386.   printf ("...sorting %d domains by area...\n", nd);
  387. #if defined(LINUX)
  388.   qsort (area, (size_t) nd, sizeof (unsigned int), (__compar_fn_t) compare);
  389. #else
  390.   qsort (area, (size_t) nd, sizeof (unsigned int), compare);
  391. #endif
  392.  
  393.  
  394. #ifdef S_DEBUG
  395.   printf ("...domain area array after sorting...\n");
  396.   for (id = 0; id < nd; id++) {
  397.     printf ("  %6u", *(area + id));
  398.     if ((id + 1) % 8 == 0)
  399.       printf ("\n");
  400.   }
  401. #endif
  402.  
  403.  
  404. /*
  405.  * for generality's sake, convert to float
  406.  */
  407.   if ((farea = (float *) calloc (nd, sizeof (float))) == NULL)
  408.       fail_alloc ("farea", 1);
  409.  
  410.   for (id = 0; id < nd; id++)
  411.     *(farea + id) = (float) (*(area + id));
  412.  
  413.  
  414. /*
  415.  * for domain coarsening data
  416.  */
  417.   if (CDYN == 1) {
  418.     for (id = 0; id < nd; id++) {
  419.       *(farea + id) /= (float) (mean_a);
  420.       if (X_SQ == ON) {
  421.         *(farea + id) = (float) sqrt (*(farea + id));
  422.       }
  423.     }
  424.     BIN_METHOD = FIX_INTERVAL;
  425.   }
  426.   ah = eval_hist (BIN_METHOD, farea, nd, 1);
  427.  
  428.   ah->mean = mean_a;
  429.   ah->sdev = sdev_a;
  430.   ah->skew = skew_a;
  431.   printf ("\n\n...histogram of %d domain areas (norm by raw mean):\n", nd);
  432.   printf ("    max(pix): %6.2f", ah->amax);
  433.   printf ("    min(pix): %6.2f", ah->amin);
  434.   printf ("      bin_w(pix): %6.3f\n", ah->bw);
  435.   printf ("  hist mean(pix): %6.2lf\n", ah->hmean);
  436.  
  437.  
  438. /*
  439.  * write histogram data to file
  440.  */
  441.   if (WRITE_FILE == ON) {
  442.     c = 0;
  443.     printf ("\n...write files (.hdt, .vin and .gdt)?(y/n)");
  444.     while ((c != 'y') && (c != 'n'))
  445.       c = readlin (in_buf);
  446.  
  447.     if (c == 'y') {
  448.       printf ("...enter fn [drive:][\\directory\\]: ");
  449.       readlin (in_buf);
  450.       strcpy (hbuf, in_buf);
  451.       strcpy (vbuf, in_buf);
  452.       strcpy (gbuf, in_buf);
  453.       strcat (hbuf, ".hdt");
  454.       strcat (vbuf, ".vin");
  455.       strcat (gbuf, ".gdt");
  456.  
  457.       if ((wfp = fopen (hbuf, "w")) == NULL) {
  458.         fprintf (stderr, "could not open output file %s\n", hbuf);
  459.         exit (1);
  460.       }
  461.  
  462.       printf ("writing data file: %s\n", hbuf);
  463.       write_aoi_data (wfp, (float) foo, nparms, ah, X_SQ);
  464.  
  465.       if ((gfp = fopen (gbuf, "w")) == NULL) {
  466.         fprintf (stderr, "could not open output file %s\n", gbuf);
  467.         exit (1);
  468.       }
  469.  
  470. #if defined(LINUX)
  471.       qsort (cbub, nd, sizeof (struct bubble), (__compar_fn_t) bcomp);
  472. #else
  473.       qsort (cbub, nd, sizeof (struct bubble), bcomp);
  474. #endif
  475.       printf ("writing data file: %s\n", gbuf);
  476.       write_bub_data (gfp, (float) foo, nparms, &cbub[0], nd);
  477.  
  478.       if ((vfp = fopen (vbuf, "w")) == NULL) {
  479.         fprintf (stderr, "could not open output file %s\n", vbuf);
  480.         exit (1);
  481.       }
  482.  
  483.       printf ("writing data file: %s\n", vbuf);
  484.       write_vin_file (vfp, nd, x1, y1, x2, y2, cbub);
  485.  
  486.       fclose (wfp);
  487.       fclose (vfp);
  488.       fclose (gfp);
  489.     }
  490.   }
  491.  
  492.  
  493.   /*
  494.    * release resources
  495.    */
  496.  
  497.  
  498.  
  499.   if (LIM_AREA == ON) {
  500.     free (rarea);
  501.     free (rcbub);
  502.   }
  503.   free (ocbub);
  504.   free (oarea);
  505.   free (farea);
  506.   free (ah->ph);
  507.   free (ah);                    /* alloc in aoi_hist */
  508.  
  509. }
  510.  
  511. /*
  512.  * extract area values into separate array
  513.  */
  514. void
  515. extract_area (int nd, unsigned int *area, struct bubble *cbub)
  516. {
  517.   int id;
  518.  
  519.   for (id = 0; id < nd; id++)
  520.     *(area + id) = (cbub + id)->area;
  521. }
  522.